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1— H , The inhomogeneous Zerilli equation is solved in time-domain numerically with the Chebyshev 

3 ■ spectral collocation method to investigate a radial-infall of the point particle towards a Schwarzschild 

' ' black hole. Singular source terms due to the point particle appear in the equation in the form of 

CNJ , the Dirac (5-function and its derivative. For the approximation of singular source terms, we use the 

direct derivative projection method proposed in [9] without any regularization. The gravitational 

waveforms are evaluated as a function of time. We compare the results of the spectral collocation 

method with those of the explicit second-order central-difference method. The numerical results 

show that the spectral collocation approximation with the direct projection method is accurate and 

^^• converges rapidly when compared with the finite-difference method. 
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^ ! I. INTRODUCTION 

J>^' The capture of a stellar mass compact object (a black hole or star) by a supermassive black hole occurs commonly 

r"| ' in the nuclei of most galaxies. As the compact object spirals into the black hole, it emits gravitational radiation which 
Qh would be detectable by space borne detectors such as NASA/ESA's LISA [ij. 

Because of the extreme mass-ratio, the small companion of a supermassive black hole can be modeled as a point 
fvj , particle, and the problem can be addressed using black hole perturbation theory. Moreover, as a first approximation, 
J> ' we will let the point particle follow geodesies in the space-time of the central black hole. This type of approach has 
l/~) ] been taken by several researchers in the past using a frequency-domain decomposition of the master equation (for 
■^ ■ a recent review see Q). That approach has had many remarkable achievements, specifically the accurate (to 10"'') 
^'J ' determination of the energy flux of gravitational waves for many classes of geodesic orbits. However, for orbits of 
^ . i: high-eccentricity and for non-periodic orbits (such as parabolic orbits, infalling or inspiraling orbits) frequency-domain 
^-H ■ computations typically result in long computation times and large errors. It is here that a time-domain based approach 
^"^ , would be invaluable. Finite-differencing based, time-domain methods for such evolutions were only recently developed 
^_,^ ' and they suffer with a relatively poor accuracy of about 1% [3, lla|- The main source of error in these methods is 
the approximate modeling of a point particle, i.e. a Dirac J-function, on a finite-difference numerical grid. Various 
approaches to this issue have been attempted, including "regularizing" the Dirac (5-function using a narrow Gaussian 
distribution [2| and also using more advanced discrete-^ models [J, [l5[ with only limited success. 

An alternate approach to the time-domain based solution of this extreme mass-ratio inspiral (EMRI) problem is to 
C^ use the spectral method instead of the finite-differencing approach discussed above. It is well known that the spectral 
method yields so-called "spectral convergence" for the approximation of smooth problems. That is, the accuracy of 
the spectral method improves by an exponential order rather than an algebraic order which is the typical order for 
the finite difference method. The spectral method, however, loses its spectral accuracy when discontinuous or singular 
problems are considered. Convergence of the spectral approximation for discontinuous problems is only 0(1) in its 
Loo norm i.e. error is 0(1) for large N if the discontinuity is included and 0(l/iV) if it is not. For the approximation 
of singular source terms, the spectral Galerkin approach yields an efficient approach because the Galerkin projection 
of the Dirac (5-function to any polynomial spaces can be easily and exactly evaluated. However, the spectral Galerkin 
approximation is only limited to 0(1) convergence near the singularity or the local jump discontinuity. The spectral 



o 



X 



'Electronic address: jaehun@buffalo.edu 
t Electronic address: ^gkhanna@uniassd.edu] 



collocation method also suffers from this 0(1) convergence and this is well known as the Gibbs phenomenon. The 
Gibbs oscillations found in spectral approximations can also easily induce numerical instability especially for hyperbolic 
equations. A regularization of the solution or the singular source terms is therefore necessary to enhance accuracy 
away from the singular regions and also to maintain numerical stability. In 9', 10] it has been shown that the direct 
collocation projection of the (5-function can yield spectral accuracy for some linear PDEs even when no regularization 
is applied. The spectral accuracy obtained by this method is due to the consistent formulation of the given PDEs in 
the collocation sense. The same direct projection method for the spectral Galerkin method, however, does not show 
such good results. Although the direct collocation method is limited only to certain classes of PDEs, it can be applied 
to more general problems if the solutions exhibit some "nice" properties such as symmetry, steady-state behavior, 
linearity, etc. @. The problem that we consider in this paper is described by the second-order harmonic equations 
with singular source terms and we will show that the direct projection method yields accurate results for this case. 

In this work, we perform a proof-of- concept computation to demonstrate that the spectral collocation method 
is more accurate and exhibits higher convergence rate when compared to a time-domain solution method for the 
EMRI problem based on the second-order finite-differencing. We demonstrate this using a non-periodic radial infall 
problem for a point particle into a Schwarzschild (non-rotating) black hole. The master equation to be solved is the 
inhomogeneous Zerilli equation [19] which is essentially a wave-equation with a potential and a complicated source 
term. We solve this equation two different ways - using standard explicit second order finite difference method and 
using a Chebyshev spectral collocation method - and compare the results. 

This paper is organized as follows. In Section 2, we briefly outline the mathematical formulation (Zerilli equation) 
behind the astrophysical process under consideration. The mass-ratio of the system is assumed to be infinitesimally 
small so that the small object can be considered as a point source. In Section 3, we explain the numerical methodology 
used for the approximation of the Zerilli equation. The Chebyshev spectral collocation method is explained in detail 
there. For the approximation of singular source terms that appear in the Zerilli equation, we introduce the direct 
collocation projection method of the ^-function and its derivative(s). Numerical comparison between the Chebyshev 
collocation method with the direct projection method and the Lax-Wendroff method for the approximation of a simple 
first-order wave equation is provided. In Section 4, numerical results of the Chebyshev spectral method are provided 
for the Zerilli equation. We also present the second-order finite difference approximations for comparison with the 
results of the spectral method. In Section 5, we summarize and suggest possible future work. 

II. EMRI USING BLACK HOLE PERTURBATION THEORY 

In this section, we summarize the equations corresponding to an extreme-mass-ratio binary black hole system. 
We will keep the discussion here very brief and simply refer the reader to the more detailed treatment presented in 
reference [1^ ■ Depicted in Figure [1] is schematic that lays out the basic set-up of the problem that we are addressing 
in this paper. The point particle initially starts out at a location tq and begins to "fall" toward the black hole event 
horizon as shown on the left side of the figure. As the particle falls, it gains more speed and radiates gravitational 
waves. Our numerical approach towards this problem uses the r* coordinate and uses appropriately chosen, finite 
values for the location of the boundaries. 

EMRI Schematic 
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FIG. 1: The schematic illustration of the EMRI for the ID inhomogeneous Zerilli equation Eq. ([T]). 7?* and ro denote the 
location of the event horizon and the initial location of the infalling point source, t denotes the physical time. As t increases, 
the infalling particle illustrated by the black dot with the arrow approaches the event horizon and emits the gravitational wave 
propagating toward _R^ illustrated by the solid ripples in the figure. 



The central Schwarzschild black hole has a mass of AI and the infalling object has a mass of uiq. For simplicity, 
the mass-ratio between the two is taken to be infinitesimally small so that the small object is considered as a point 
source moving on a geodesic of the curved space-time due to the other, i.e. uiq/M —> 0. Because we are treating 



this problem in the context of first-order perturbation theory, the governing equation of this system is hnear and 
the solution is the same for any mass-ratio. One solution of a particular mass-ratio can be derived from another 
of different mass-ratio, simply by rescaling. We assume that the point source is falling along the z-axis toward the 
hole without any angular momentum. As the point source is falling into the massive hole, gravitational waves are 
emitted by the system. The lowest order perturbation theory of the initial Schwarzschild black hole spacetime leads 
to the inhomogeneous Zerilli equation with even-parity. Such an equation describes the gravitational wave i/) in 1 + 1 
dimension given by the following second-order wave equation with a potential and source term, 

i^u = i^r'r'-Vi{r)^-Si{r,t), (1) 

where r* is the tortoise coordinate, Vi{r) the potential term and Si{r,t) the source term [i2,lly]- Note that we have 
set the units c = G = 1. The tortoise coordinate mentioned above, is given by 

r* ^r + M\TL{r/2M ~l), (2) 

where r is the standard Schwarzschild coordinate and r > 2M . Note that tortoise coordinate r* approaches —oo as 
the value of r approaches the Schwarzschild radius r = 2M. The Zerilli potential Vi{r) is given by 

^^^'^ = f^'^ .3(Ar + 3Af)^ ' (') 



where the function /(r) is 



and the parameter A is 



/(.•) ^ 1 - ^ (4) 

r 



X = {l + 2){l^l)/2. (5) 



In our case of interest, Z-pole is only even with azimuthal angular term m = 0. It is noteworthy that as r —^ 2M , f{r) 
becomes /(r) -^ and the potential vanishes. As r ^ oo, the potential also vanishes. The local maximum of Vi{r) 
exists near the Schwarzschild radius. 
Now, the source term is given by 



2{l-2M/r)K 
Si(r,t) = - 



-r^l^2M/r)^S'ir-R) 



r{X + l){Xr + 3M) 

r{X + 1) - M 3Af [/°r(l - 2M/rf 



2[/o rX + MI 

where k is 

K = 87rTOo\/(2^ + l)/47r, 
and the zero component of the 4- velocity of the particle C/" is given by 



6ir - R) 



(6) 



0^ / l-2M/ro 
Y 1 - 2M/r ' 

where tq is the initial position of the particle. 

Singular source terms S{r — R) and 5'{r — R) are due to the point source and its movement along r. Here R is 
the coordinate of the point source at time t. Although the source term expression takes a complicated form, one can 
easily obtain some intuitive aspects of the physical system from it. For example, note that the term simply scales with 
mo (through the quantity defined above as k) which is simply a reflection of our linear approximation as mentioned 
earlier. In addition, as expected, the expression is non-zero only at the location of the point particle (r — R). It 
is also interesting to note that, if the particle is very far from the black hole, it is the 6'{r — R) term that strongly 
dominates the source expression, while in the other extreme case the entire term becomes negligible. The superscript 
' denotes the derivative with respect to r not r*. Since the derivative with respect to space in Eq. ([1]) is given in r* 
rather than r, we need to convert r to r* to locate the point source at R. For our numerical approximation in the 



next section with the spectral method, r* is transformed to the unit interval [—1,1] using the linear map as explained 
in the next section. The free fall trajectory R{t) of the test particle is then given by 



If To — + cx), the initial system can be regarded as the unperturbed Schwarzschild spacetime and the wave function ip 
vanishes for all r at t = 0. Due to the computational limit, however, ro is only finite for the numerical approximation. 
For the current work, we assume that mo is infinitesimally small compared to the mass M of the Schwarzschild 
black hole so that the initial ip is chosen to vanish. This causes a small discrepancy in the full solution. The initial 
fluctuations due to such discrepancy, however, leave the computational domain according to Eq. ([1]) as t — > oo. 

III. NUMERICAL METHODOLOGY 

In this section, we present the numerical methodology used in this paper - the spectral collocation method and the 
direct projection method of the (5- function to solve the second-order wave equation given in Eq. ([1]). A numerical 
example for the simple advection equation is also provided. For the time integration, the TVD third-order Runge- 
Kutta method is used for the spectral method. 

A. Chebyshev spectral collocation method 

For the spatial differentiation, the Chebyshev spectral collocation method is used. Suppose that the differential 
operator in one dimension is given on the solution u{x, t) for x £ Q and t > as 

Cu = 0, (8) 

with the boundary conditions 

B^u = h^{t), xedn (9) 

where B^ are the boundary operators at the domain boundary dfl. Let Vn be the projection operator applied to u{x, t) 
onto the polynomial space Bn expanded by the Chebyshev polynomials of up to degree A^. The Chebyshev spectral 
collocation method seeks a solution in the Chebyshev polynomial space expanded by the Chebyshev polynomials T/ (x) 
as 

m 

UN = rNu{x,t) = Y,Mt)Ti(x), (10) 

where Ti(x) is the Chebyshev polynomial of degree I and ui the corresponding expansion coefficient. Chebyshev 
polynomials are defined by the following orthogonality property 

1 '■' 



w{x)Ti{x)Ti'{x)dx = Sw, (11) 

hi J_i 

where hi is the normalization factor given hy hi — n for I = and hi — ^ for / G Z+ and w{x) is the Chebyshev 
weight function given by w{x) = 1/Vl — x^. The expansion coefficients {ui}^^ are found based on approximations 
at the collocation points defined in CI. Here Sw is the Kronecker delta in a usual sense. Then the exact expansion 
coefficients ui is obtained by 

1 '■' 



ui{t)^— w{x)u{x,t)Ti{x)dx. (12) 

hi J-i 

Since the collocation method seeks the solution at the particular collocation points, ui(t) is obtained by evaluating the 
integral using the quadrature rules on such collocation points. As long as the notation is not confusing, we also use 
ui as the approximation of the exact expansion coefficients. Since the boundary conditions are given, the commonly 



used collocation points are the Gauss-Lobatto quadrature points. The Chebyshev Gauss-Lobatto collocation points 
Xj are given by 

x, --cos(^j), Vj = 0, •••,]¥, (13) 



here N — m for the completeness. With this condition, we construct the approximation based on the solution at the 



collocation points, {u{xj^ty\^^Q such as 



*(a;, ^) = X! ^'(*)^Ka;) = X! V'Ka;)M(2;i, = X! '^'(^)"'' (1^) 



i=0 1=0 1=0 



where let ui denote the grid function of u{x,t) at a; = Xj and ipiix) are indeed interpolating polynomials. The 
expansion coefficients are given by 



2 ^ 1 

^'W = ^^Er^^w^'(^^)- (15) 

' i=o -J 

The interpolating polynomial ipi{x) is a (^-function in the sense that ipi{xi') = Sw. One may adopt 4'i{^) &s an 
approximation of the (5- function, but using ipii^) as a (5- function for the given PDE yields only highly oscillatory 
solution near the singularity [9]. The interpolating polynomials ipii^) on the Gauss-Lobatto-Chebyshev collocation 
points are given by 

where Q = 1 for I = 1,- ■ ■ , N — 1 and ci — 2 for I ~ 0,N. By taking a derivative of ipi{x) we obtain the derivative of 
u{x,t) with respect to x such as 

rn 

u'(x,t)=^^;(x)uj, (17) 

3=0 

where the superscript ' denotes the derivative with respect to x. Then at the collocation points we obtain 

u' = D • u, (18) 

where u = (uo, • • • , un)'^ and D the differentiation matrix given by the derivative of the interpolating polynomial at 
the collocation points as 

Di,^^j'^{xi). (19) 

Thus we know that every column of D is the derivative of the (5- function at the collocation points. The elements of 
D with the Chebyshev Gauss-Lobatto collocation points are given by [0,3, 

A, = <! ^1-? ' = '■' '^°'^ . (20) 

I — 

i^ N 

The Chebyshev collocation method seeks a solution u(x,t) by having the following discrete residue Rn{x) vanish 
at the collocation points, 

Rn{x) ^ Cnun -0, X ^ Xj,yj ^0,- ■ ■,N, (21) 

where Cn — Vn^Vn- 
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B. Spectral projection of singular source terms 

For singular source terms, we use the direct projection method 9]. The direct projection method is to approximate 
the ^-function in a consistent way with the given PDE. Suppose that we are given a PDE in one dimension 

where C = d/dx and the boundary condition u(— 1) = 0. Here the (5-function in the bounded domain Q is defined as 

S{x)dx = 1, and / S{x)u{x)dx = ^(0). (22) 

-1 J-i 

Then the exact solution is given by the Heaviside function Hc{x) 

where c — for this case. Then the direct projection method for the (5-function is obtained by simply mimicking the 
relation 

A.H,{x)^6{x), 

at the collocation points such as 

dN = B U, (23) 

where 5n is the collocation approximation at the collocation points, Sn — {Sn{xo): ' ' ' iSn{xn))'^ , and H is H = 
(0, • • • , 0, 1, • • • , 1)"^. For H the nonzero value begins at x = Xj such that Xj > c. The consistent formula of 6^ does 
not necessarily satisfy the definition of the (5-function in Eq. (j22p . Suppose that u{x) = |s(7n(a;)|. Then by definition 
at each collocation point Xj , 

d r 

— / S{s)u{s)ds = S{x)u{x) = 0, for Xj,\fj — 0,- • • , N. 

At Xj the consistent formula of 5n yields 

5{x)u{x)\^^ ^ (D • IiN)ju{xj) = (D • H)j ^ 0, for a;^, Vj = 0, • • • , iV. 

For the time dependent problem such as the Zerilli equation considered in the next section, the source term is also a 
function of time. In this case, the parameter c for Hc{x) is a function of time, i.e. c = c(i). If c{t) coincides with one 
of the collocation points, Xj at a certain time tp, the definition of the Heaviside function Hc{x) at Xj and tp can be 
either iJc(t ){xj) = 1 or /^^(i ){xj) = \- Both definitions yield the same results. 

C. Numerical examples 

To demonstrate the accuracy obtained by the direct projection method for the ID wave equation with the singular 
source term, we consider the following advection equation 

ut + Ux = 5{x), a;e[-l,l], f > 0, (24) 

with the initial condition u{x^Q) — — sin(7ra;) and m(— l,t) = — sin(7r(— 1 — t)) for i > 0. The exact solution of Eq. 
([24ll is then given by 

u(a;,i) = — sin(7r(a; — t)) + / 5{x)dx. 

As the exact solution implies, the solution u{x^ t) is discontinuous at a; = 0. The jump magnitude [u] at a; = is 1 for 
any time t, 

MU=o := m(0+) - u(0") = / 5{x)dx - I 5{x)dx = 1. 



and it is clear that 



For the numerical experiment, we consider N = odd. Then the exact solution u{x, t) at the Gauss-Lobatto collocation 
points X = Xi and at t are given by 

r -sinWx.-i)) zg[0,...,i^] 
^^'' ^"11- sin(^(x, - t)) * G [^, • • • , TV] ' 

with a;(Ar_|_i)/2 < < a;(jv_|_3)/2- Let v be the spectral approximation to Eq. ([M)) such that u — (wq, • • • , v^)^ ■ Here Wj 
denotes the grid function at a; = Xj. Then we seek v with the direct projection method such that 

on 

— + D-w = D-H, (25) 

dt 

where ^ denotes the time integration operator applied to the discrete v. For the time integration, we use the 3rd 
order TVD Rungc-Kutta scheme [lj| . The boundary condition is directly applied 

Wo = -sin(7r(a;o - i)), (26) 

after the final stage of the RK step at the given time t. 

Suppose that the numerical approximation of v(x, t) is given in the form of v{x, t) = w{x — t) + z{x) where w{x — t) 
and z{x) are the time-dependent and time-independent parts of v{x,t) respectively. With the method of line, Eq. 
becomes 

^ + D • (z« + z) = D • H, 

dw .^ „ .... ..-. ^., 

hDw = 0, Dz = DH. 

dt 

The boundary condition is given by 

v{-l, t) = u(-l, <) = - sin(7r(-l - t)) + / 5{x)dx ^ - sin(7r(-l - t)), 

that is, w{—l, t) = — sin(7r(— 1— i)) and z{—l) = 0. Thus we know that 'w{x^ t) is the spectral collocation approximation 
oi wt + Wx =0 because its exact solution is given by w{x, t) — — sin(7r(x — t)). The time-independent solution z{x) is 
then obtained from by solving D • z = D • H. Since the boundary value of z{x) is given as z{—l) — 0, z{x) is given by 

D- (z-iJ) = 0, 

where the superscript " denotes the inner matrix or vector without the first column and row. D is singular, but D is 
non-singular. Thus z{x) = H{x) exactly. Here note that such exactness is obtained only at the collocation points. 
Due to the uniqueness of the solution of Eq. (|25p . v{x, t) = w{x, t) + z{x) indeed. The error function E{x, t) is given 
by E(x, t) = u{x,t) — v(x,t) = — sin(7r(x — t)) + H{x) — {w{x,t) + z(x)) = — sin(7r(a; — t)) — w{x, t). Thus we know 
that the error is solely determined by the homogeneous solution w(a:, t) which is smooth function for Vt > 0. Thus, 
the error decays spectrally. 

Since the initial condition v{x,Q) is simply given by w(a;,0) — — sin(7r(a; — i)), there exist the Gibbs oscillations 
due to the initial discrepancy of the initial condition. To see how the initial discrepancy changes with time, i > 0, 
let e(a:,t) = v{x,t) — v'^{x,t) where v'^(x,t) is the numerical solution with the exact initial condition. By definition, 
e{x, 0) 7^ 0. Since the boundary condition is the same for both v and v'^, we know that e(— 1, t) — for Vf > 0. Since 
both V and v'^ are obtained from the same equation, we obtain 

de „ 

Since e(— l,t) = for Vi > 0, e{x,t) = for t > 2. Thus the initial Gibbs oscillations due to the discrepancy of the 
initial condition do not affect the solution behavior after a finite time t. 

Finite difference scheme: For the comparison, we consider the Lax-Wendroff scheme to solve Eq. (|24[) with the 
direct projection method. Suppose that the finite difference solution is sought with the evenly spaced grids for a given 
A^ such that 

Xj = xq + jh, xq — —1, h = Ax — 2/N. 



Let D~^ , D~ and D^ the forward, backward and central difference operators defined as [31 

hD^Vj = fj+i — Vj, hD^Vj = Vj — Vj-i, 2hD Vj = fj+i — I'j-i- 

We also use the direct projection method for the approximation of the (5-function as 

dHN(x) 

On{x) = — , 

ax 

where -^ is now the central difference operator D^ with the Lax-Wendroff scheme such that 

dHMJx) Hm{Xj^i)~Hn{Xj-i) 

-^—-.DH^ix)^ , (27) 

where H]^(x) is the grid function with the finite differencing grid of Hq{x). This definition of S]^(x) satisfies the 
definition of the (5-function almost everywhere in the domain. For u{x) = \sgn{x)\, Sn{x) yields for even N, 

Then the Lax-Wendroff scheme is given by 

7,2 



where dt is the time step and we choose dt as small as the stability is guaranteed. 
To measure the errors we define the discrete L2 and Loo errors as 



L2 Error 
and 



1 ^ 

\ 1=0 



Loo Error — max \v{xi,t) — u(xi,i)|. 

ie[0,---,7V] 

Figure [H shows the approximations Sn of S{x) with the spectral method given in Eq. (|23p with iV = 65 (left) and 
the centered difference scheme given in Eq. (l27l) with iV = 66 (right). As shown in the figures, the direct projection 
approximation of 6{x) with the Chebyshev method is highly oscillatory over the whole domain. Note again that N is 
assumed to be odd to avoid any ambiguity at a; = c = for the spectral method. The approximation with the centered 
difference method, however, yields the approximation without any oscillations unlike the spectral approximation. The 
figure also shows that the J-function is highly localized around the singularity at x = with the finite difference 
method. The oscillations shown in the spectral approximation is simply, the Gibbs phenomenon. It is also well 
known that the spectral solution is highly oscillatory if there is a local jump discontinuity such as the (5-function. As 
mentioned earlier, in Q it has been shown that the direct projection method yields an accurate solution and spectral 
accuracy can be also achieved for certain types of PDEs although the approximation of the (5-function obtained by 
this method is highly oscillatory. 

Figure [3] shows the solutions and pointwise errors a.t t — 10. The left figure of Figure [3] shows the solutions with the 
Chebyshev spectral (top) and Lax-Wendroff (bottom) methods. The red and blue lines represent the exact solution 
and the approximation, respectively. For the spectral approximation, the exact solution and the approximation are 
not distinguishable. There are no Gibbs oscillations observed in the spectral approximation. It seems that the direct 
spectral approximation of the (5-function yields a cancellation of these oscillations. In [9l|, it has been discussed that 
such cancellation is due to the consistent formulation to the given PDF. For the approximation with the Lax-Wendroff 
method, small oscillations are seen around the singularity at a; = 0. In the right figure of Figure[31 the pointwise errors 
versus x are given. The figure shows that the Chebyshev spectral method yields very accurate result over the whole 
domain without any indication of the degradation of convergence order. Indeed, with the Chebyshev spectral method, 
the spectral accuracy and uniform convergence are recovered even with the local jump discontinuity. Contrary to the 
spectral approximation, the approximation with the Lax-Wendroff method only shows the low order accuracy. 

Figure 0] shows the trace of approximations with time with the Chebyshev spectral method (left) and the Lax- 
Wendroff method (right). There are small oscillations in the beginning time stages with both methods. These 
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FIG. 2: Left: The direct spectral projection approximation of S{x) on the Chebyshev Gauss-Lobatto collocation points. Right: 
The direct central difference approximation of 5{x) on the evenly spaced grids with the central difference operator D" forN = 66. 




- 



fHrfp^y~'''''T^^ 



2 , Lax-Wendroff 




FIG. 3: Left: Solutions att = lQ with the Chebyshev spectral (top) and Lax-Wendroff (bottom) methods. The red and blue lines 
represent the exact and the numerical solution, respectively. Note that the exact and numerical solutions are not distinguishable 
for the Chebyshev spectral method. Right: Pointwise errors Log\o\u{x,t) — v{x,t)\ versus x with the Lax-Wendroff method and 
Chebyshev spectral methods. 

oscillations are due to the inconsistent initial condition. These oscillations with the Chebyshev spectral method, 
however, leave the domain quickly as time goes on while they remain with the Lax-WendrofF method. The right figure 
for the Lax-Wendroff method clearly shows that such oscillations exist with all time t S (0, 5] while no significant 
oscillations are seen in the left figure for the Chebyshev spectral method even after t > T. Here T is the time for the 
initial Gibbs oscillations to leave the domain and it is about T ^ 1 as the wave speed is 1. 

Table 1 shows the L2 and Loo errors of the solution for various approximations of the (5-function. In the table, LW- 
Direct, LW-Gaussian, SP-Direct and SP-Gaussian represent the Lax-Wendroff method with the direct and Gaussian 
approximations of the 5-function, and the Chebyshev spectral method with the direct and Gaussian approximations 
of the (5-function, respectively. As shown in the table, the best result is obtained with the Chebyshev spectral method 
with the direct approximation of the (5-function. The table also shows that the direct projection method yields a 
slight better result for the Lax-Wendroff method. 



IV. NUMERICAL APPROXIMATION OF ZERILLI EQUATION 



Now we consider numerical approximations of the inhomogeneous Zerilli equation with the spectral and finite 
difference methods. For both the spectral and finite difference methods, we use the direct projection method of the 
5-function. Here we note that the finite difference method used in this work is the explicit second-order centered finite 
difference method. High order finite difference methods for the Zerilli equations such as a fourth-order convergent finite 
difference method, have been developed |ll| . The comparisons between the spectral and high order finite difference 
approximations for the Zerilli equations will be left for our future work. 
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FIG. 4: Trace of solutions with time with the Chebyshev spectral (left) and Lax-Wendroff methods (right). 



TABLE I: Z/2 and Loo errors with various approximations of 5{x) at t = 10 for Eq. 



Method 



1/2 error Loo error 



LW-Direct 1.02(-2) 2.52(-2) 

LW-Gaussian 4.52(-2) 2.58(-l) 

SP-Direct 3.60(-10) 7.77(-10) 

SP-Gaussian 4.04(-2) 2.33(-l) 



(n) = 10" 



The governing equation is defined on the simple 1 + 1 spacetime with time and space denoted by t and r* respectively. 
The range of r* , the tortoise coordinate is given by 

K <r*< R*^, 

where i?* is obtained close to the radius of the event horizon and i?^ — > oo. In practice, i?* and R^ can not be — oo 
and oo respectively. Let i?* is taken at r ~ 2M and R%^ at r 3> 2M. By definition, i?* < and i?^ > 0. The tortoise 
coordinate r* is linearly mapped into the Chebyshev Gauss-Lobatto collocation points {£,j}jLo given as 



where 






C,=-cos{7TJ/N), j = 0,---,N. 



A. Derivative operator D 



(28) 



Using the linear mapping of Eq. ([28]) , the derivative operator with respect tor*, D is given by 



D:-- 



d 2d 



dr* Rl, - Rl 9e Rio - R. 



-DiO, 



where D{^) is the differential operator defined on the Gauss-Lobatto collocation points {^j}f^Q given in Eq. (|20p . 
The second derivative operator is given by 

Then the Chebyshev spectral method for the black hole system Eq. ^ seeks an approximation for 



dt^ 



D'i; - Vi^P - S, 



(29) 



11 

where -0, Vj and 5* are the wave function, the potential of I mode and the source term defined on {CjJwLoj respectively. 
Singular source terms S{r,t) and S'(r,t) are approximated with the direct projection method as 

S{r-R)^D- Hr, 6'{r - R) ^ D ■ 5{r - R) ^ D^ ■ Hr, 

where H is the Heaviside function defined on {£,j}f^Q- Here we note that i? is a function function of time and the 
Heaviside function Hr is being shifted as time goes on. 

In our computation, the second derivative of ip is with respect to r* not r. Thus at each collocation point, r*, the 
corresponding proper distance should be calculated to obtain the potential and singular source terms accordingly, for 
which the bi-section method is used. 

B. Time stepping 

For the time-stepping, we adopt the second-order time stepping method 

i;{t + At) - 2-ip{t) + ipU - At) 
^t^ = A^2 • 

This method is two-step method. Thus we need two initial conditions for the time marching. We can simply use 

V^(r,0)=V(^Ai), 

or we can evolve Eq. (1) using the first order approximation for ip{r, At) [7|. 

C. Initial condition 

The initial location of the point particle is r = ro < oo. We have to find the exact initial condition when the particle 
is located at finite r — rg initially. In this paper, we assume that 

V'(r,t = 0) = 0. 

This vanishing initial condition is not physical and inconsistent to Eq. ([1]) as it does not satisfy the governing equation. 
Due to the inconsistent treatment of the initial condition, the numerical approximations of ^ in the beginning stages 
are highly oscillatory and it takes a certain time for these initial oscillations to leave the computational domain as we 
will see in the next sections. 

D. Boundary conditions 

We impose the absorbing boundary condition at the boundaries r* = R* and r* = R*^. Since this is the ID 
problem, the following boundary condition should be a perfectly absorbing boundary condition: 

V^t-V-r-^O, at r*=Rl, 

and 

V't+^r'-O, at r*=R*^. 

Here we note that the derivative with the spectral method at the end points should not be obtained using the global 
derivative operator D because there can be a discontinuity inside the domain. That is, the first order derivative term 
is not given as D ■ ip\r=R^- The no-flux boundary conditions should be only local. We use the upwind method at both 
ends such as 

Hr*o,t + At)-^irlt) i,{rt,t)-i,{r*o,t) , 

TT = ; ; > at r = R^, 

At r^ - To 

where Tq = i?* and 

V'(rS,t-f-At)-7/;(rc5,t) i^{r*f„t) - i:ir*^_^,t) 



At r% - r^_i 



at r* = R*^, 



where Tq = i?* 
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E. Spectral filtering methods 

We also adopt the spectral filtering method to minimize the possible non-physical high-frequency modes. The 
oscillations with the spectral method possibly found near the local jump discontinuity and also generated due to the 
inconsistent initial conditions, propagate through the whole domain. To reduce these non-physical modes, a filter 
function is applied to the solution. To apply the filter function, we transform the approximation found in the physical 
domain to the coefficient domain using Eq. (jlSp . Then we obtain the filtered approximation as 



N 



3=0 



where ail/N) is the filter function according to 20]. The filtering method is equivalent to the spectral viscosity 
method [8|. For the current work, we use the exponential filter defined as 

<j{j/N) = exp(ln eij/N)P), Vj = 0, • • • , iV, 

where e is a strictly positive real constant and p is known as the filtering order. The filter function is constructed to 
satisfy the properties, cr(0) = 1, o'^''^(0) = for g = 1, • • • ,p — 1, cr(±l) = 0, and a{—ri) = a{r]) [23|- In principle, 
ujf{t) -^ UN{t) if p ^ oo. If p ^ 0+, ujf{t) becomes only 0th order approximation of UAr(i). By the definition of the 
filter function (j{j/N), cr(l) -^ 0, for which the positive constant e is usually chosen so that exp(lne) becomes about 
machine accuracy, eM- For our work, we choose e^v/ ~ 10^^^. Using Eq. (|15p . we have 

N 

u% {xi) = ^ SijUj , Vi = 0, • • • , iV, 
j=o 

where Stj are 

1 ^ 2 
^''^- = - E ^Tn{x,)T^ixj) exp(lne(n/iV)P). 

•' n—0 

F. Filtering for the finite difference method 

For the finite difference method, we also experimented with a filtering method. We simply convolved a Gaussian 
distribution with the source-term itself, yielding a smoother profile for it. This Gaussian filter approach has been 
used elsewhere [2l| in literature is a natural one to attempt in our context. Through experimentation, we realized 
that choosing the width of the Gaussian to be of the same scale as the spacing between grid points yields the best 
results. Figure [7] depict sample results from performing such Gaussian filtering in the finite difference evolutions. 
These results show that the Gaussian filtering yield better results than those without the filtering. With the Gaussian 
convolution the filtered function fj^ '^^'^ for /jv is given by 



j:fxitered^^^_l_-^ ^ / f Niv ,t) exp{-{x - yf / a)dy , 



IN 

The Gaussian filter is applied for the source term. Since the source is only localized, the integration is done only over 
a certain number of grids, that is, for example at ith grid point, x = Xi, the integration on the grids is done over 
[xi-n , Xi+n ]• Thus totally 2np -I- 1 grid points are used for the integration. According to our numerical experiments, 
the best parameters for the Gaussian kernel for various N are {N,np,a) — (32, 2,8), (64,6, 8), (128, 32, 16), and 
(1024,32,32). These values are used for the filtering for the finite difference method in the following sections. 

G. Numerical results 

For the numerical experiments, we use M = 1, mg = lO^^*', i?* — —30, i?^ = 30, 130, and I — 2. For the time 
integration we use dt = 5x 10~^ with the final time Tfinai — 300, 600, and 2400. For the spectral case, filtering orders 
in the range of p = [32, 64] were used. For the finite difference method, the filtering parameters given in Section 4.6 
are used. For the comparison of the spectral method with the finite difference method, we focused on accuracy and 
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convergence of the ringdown profile with time and N. Since the current work is based on the spectral method with the 
single domain, the finite difference method has faster computing time than the spectral method. The multi-domain 
spectral method will be considered in our future work to increase the computing speed. 

Figure [5] shows the waveform ip/rriQ of unfiltered and filtered spectral solutions for 130 < t < 180 with N = 128. 
For filtered solutions, p — 64,48, and 32 are used. The left figure shows the case of filter on-set time, tfuter = and 
the right of tfuter — 10. The oscillations shown in the unfiltered solution are considerably reduced in the filtered 
solutions, and the filtered solutions with different p agree well. It is also observed that the filter on-set time does not 
affect the solution in the ringdown phase if tfuter < t '^ 130. 

Figure [6] shows the waveforms ip/mo of filtered spectral solutions and those of the finite difference method. N = 128 
is used for the spectral method and N — 128 (left) and N = 1024 (right) are used for the finite difference method. 
No filtering is used for the finite difference method. The figures show that the waveforms are highly oscillatory with 
the finite difference method and also show that more oscillatory waveforms are found with larger N. Furthermore 
the figures show that the global behavior of the finite difference solution with large N = 1024 agree with that of the 
filtered spectral solutions better than with N = 128. This implies that the spectral method can yield more accurate 
solution with smaller N than the finite difference method with the proper choice of p. 





FIG. 5: xp/mo versus t with the spectral method for p = 64 (red), p = 48 (blue), p = 32 (magenta) and unfiltered method (black) 
for t G [130, 180]. Left: tfuter = 0. Right: tfuter = 10. 




Finite Difference 
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Finite Difference 
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FIG. 6: xp/mo versus t with the spectral and finite difference methods. For the spectral method, p — 64 (red), p = 48 (blue), and 
p — 32 (magenta) are used with N = 128 and tfuter ~ for both plots. For the finite difference method, N = 128 (left) and 
N = 1024 (right) are used. 



Based on the solutions obtained, we conducted a Fourier analysis to find the period T and the frequency uj of the 
ringdown waveforms. Table HIl shows the period and frequency obtained with the Fourier analysis. In the table CD 
denotes the finite difference method with the central differencing, SP the spectral method and p the spectral filtering 
orders. The period computed based on the frequency domain approach is well known as T ~ 16.815. The table clearly 
shows that the periods of the ringdown waveforms with the spectral method for TV = 32 to iV = 128 and with the 
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finite difference method for N = 128 and N = 1024 are about T ^^ 16.8. The finite difference method with smaller N 
such as N — 32,48, and 64, however, yields inaccurate results of T ~ 18.7 and T ~ 15.3. The table also shows that 
the spectral method yields an accurate period even with smaller N such as A^ = 32. 

TABLE II: The period T and frequency lo of the ringdown of the waveform ip/mo for 1 = 2. 



Method 


N 


P 


U! 


T 


CD 


32 




0.33649906435356 


18.672 


CD 


48 




0.33649906435356 


18.672 


CD 


64 




0.41127663420991 


15.277 


CD 


128 




0.37388784928174 


16.805 


CD 


1024 




0.37388784928174 


16.805 


SP 


32 




0.37388784928174 


16.805 


SP 


48 




0.37388784928174 


16.805 


SP 


64 




0.37388784928174 


16.805 


SP 


128 




0.37717970828492 


16.658 


SP 


128 


64 


0.37388784928174 


16.805 


SP 


128 


48 


0.37388784928174 


16.805 


SP 


128 


32 


0.37388784928174 


16.805 



Figure [7] shows the waveform in logarithmic scale for t e [2000, 2070]. The left figure shows the unfiltered spectral 
approximations with N = 128 (black), A^ = 64 (red), iV = 48 (blue), and TV = 32 (magenta). The right figure shows 
the filtered finite difference approximations with N — 1024 (black), N — 128 (red),iV = 64 (blue) and TV = 32 (green). 
The figures clearly show that the spectral approximation yields a fast convergent behavior as N is increased, without 
any filtering method applied. The finite difference approximations, on the other hand, yield only a slow convergence 
as the right figure indicates. Moreover, large iV such as TV = 1024 is needed for the finite difference method to achieve 
comparable accuracy as the one obtained by the unfiltered spectral solution with N = 64. It is worth noting that the 
convergence with the finite difference method is possible only if the Gaussian filtering is used. 
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FIG. 7: Logioltl^/mol versus t. t £ [2000, 2070]. Left: Unfiltered spectral approximations with N = 128(black), QA(red), 4:8(blue), 
and 32 (magenta). Right: Filtered finite difference approximations with N = 1024 (black), 128 (red), 64 (blue), and 32 (green). 



SUMMARY & OUTLOOK 



In this work, the inhomogeneous Zerilli equation has been solved numerically using both the Chebyshev spectral 
collocation method and the 2nd order explicit finite difference method to investigate the radial-infall of a point particle 
into a Schwarzschild black hole. The infalling object is considered as a point source and consequently singular source 
terms appear in the governing equations. For the approximation of the singular source terms, we use the direct 
derivative projection of the Heaviside function on the collocation points. Such an approach yields very oscillatory 
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approximations for the (5-function and its derivative. Despite these oscillations, the spectral solution is obtained in 
an accurate way as the direct derivative projection is consistently formulated with the given PDEs. A simple wave 
equation in 1 + ID spacetime shows that the spectral method recovers the exponential convergence even with the local 
jump discontinuities. The same technique is applied to solve the Zerilli equation and various numerical results have 
been presented. These results are also compared with those from the second-order explicit finite difference method. 
The spectral approximation with large A^ shows small oscillatory behavior after the effect of the initial oscillations is 
minimal. To reduce these unphysical oscillations, the spectral filtering method is applied. A proper order p of filtering, 
yields a smooth waveform with these oscillations considerably reduced. The finite difference method needs large A^ 
to obtain the approximation with the same accuracy of the spectral method. Therefore, the spectral method yields 
accurate results with much smaller A^ and can be efficiently used for the approximation of these types of problems. 

Our future work will center around the development of more efficient spectral methods such as those with the domain 
decomposition techniques and also adaptive filtering. These will further improve the accuracy and the computational 
efficiency of the solution. For the multi-domain spectral method or the discontinuous Galerkin method, it is possible 
to place the singular source term at the domain or element interface so that the singular source term plays a role 
only as a boundary or interface condition. In this approach, the singular source term is involved in the given PDEs 
only in a weak sense. For the spectral multi-domain penalty method, the penalty parameter can be also exploited to 
enhance both stability and accuracy associated with the singular source term. In general, however, the location of 
the singular source term varies with time and it would be more convenient to allow the singular source term to exist 
inside the sub-domain or element, in particular when the range of the possible locational variation is large. 

Finally, we will also investigate more general collisions of black holes (for example, in 2 or 3 dimensions, with 
rotation, etc.) and compare with the results from the frequency domain method. 
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